########## re run everything
rm(list=ls())

# load packages
library(ggplot2)
library(ggbreak) 
library(readxl)
library(dplyr)
library(tidyr)
library(stargazer)
library(texreg)
library(sjPlot)
library(estimatr)
library(sandwich)
library(lmtest)
library(miceadds)
library(lme4)
library(lmerTest)
library(glmnet)
library(plm)
options(scipen=999)


# load data
# set working directory to replication file
load("~Damann Kim Tavits/6. Submissions- IO/5. Replication/PostUnit.Rdata")
load("~Damann Kim Tavits/6. Submissions- IO/5. Replication/PoliticianDayUnit.Rdata")

# number of politicians
n_distinct(polday$id) # 469
n_distinct(post$id) # 469



### Table S2.1
topic_dist <- post %>% group_by(label) %>% summarise(prop = round(n()/nrow(post)*100, digits=1))
print(topic_dist)



# create a pdataframe for panel regressions
panel <- pdata.frame(polday, index = c("id", "date"))

# Table S3.1 (Controlling for Political Office)
mod1fe <- plm(num.post ~ daysince1 + invasion*woman + daysince2 + office, 
              data=panel,
              model="random", index = "id")
nwse <- coeftest(mod1fe, vcov.=vcovNW)[,2]
pval <- coeftest(mod1fe, vcov.=vcovNW)[,4]

screenreg(mod1fe, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Posts per politician",
          custom.coef.map = list("invasion1"="Invasion",
                                 "woman1"="Woman",
                                 "invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "daysince1"="Time",
                                 "daysince2" = "Time since invasion",
                                 "officegovernment_cabinet" = "Cabinet",
                                 "officehead_regional_state_administration" = "Governor",
                                 "officemayor" = "Mayor",
                                 "officemp" = "MP",
                                 "(Intercept)"="Intercept"),
          caption = "Effect of War on Politicians' Engagement with Office FE",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:office",
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod1fe$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469),
          reorder.gof = c(1, 3, 2))


#Table S3.2: Controlling for Regional Variations
mod1region <- plm(num.post ~ daysince1 + invasion*woman + daysince2 + region, 
                  data=panel,
                  model="random")
nwse <- coeftest(mod1region, vcov.=vcovNW)[,2]
pval <- coeftest(mod1region, vcov.=vcovNW)[,4]

screenreg(mod1region, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Posts per politician",
          custom.coef.map = list("invasion1"="Invasion",
                                 "woman1"="Woman",
                                 "invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "daysince1"="Time",
                                 "daysince2" = "Time since invasion",
                                 "(Intercept)"="Intercept"),
          caption = "Effect of War on Politicians' Engagement with Region Controls",
          label = "tab:region",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod1region$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469,
                                 "Region controls" = "\\checkmark"),
          reorder.gof = c(1, 3, 4, 2))


# Table 3.3: Distribution of Committee Assignment by Gender
dist <- panel %>% group_by(committee) %>%
  summarise(n = n_distinct(id),
            n.woman = n_distinct(id[woman==1]),
            prop = n.woman/n) %>% distinct() %>%
  arrange(desc(prop)) %>% mutate(prop = round(prop, digits=2))
print(dist, n=25) 


# Table 3.4: Committee Assignments and Posting Behavior
mod1com <- plm(num.post ~ daysince1 + invasion*woman + daysince2 + committee, data=panel,
               model="random", index = "id")
nwse <- coeftest(mod1com, vcov.=vcovNW)[,2]
pval <- coeftest(mod1com, vcov.=vcovNW)[,4]

screenreg(mod1com, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Posts per politician",
          custom.coef.map = list("invasion1"="Invasion",
                                 "woman1"="Woman",
                                 "invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "daysince1"="Time",
                                 "daysince2" = "Time since invasion",
                                 "(Intercept)"="Intercept"),
          caption = "Effect of War on Politicians' Engagement with Committee Controls",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:committee",
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod1com$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469,
                                 "Committee controls"="\\checkmark"),
          reorder.gof = c(1, 3, 4, 2))

# Table S3.5: Foreign Policy Elites and Posting Behavior
panel <- panel %>% 
  mutate(foreign = as.factor(ifelse(committee %in% 
                                      c("Ukraine's Integration into the European Union",
                                        "National Security, Defence and Intelligence",
                                        "Foreign Policy and Inter-Parliamentary Cooperation"),
                                    1, 0)))
mod1for <- plm(num.post ~ daysince1 + invasion*woman + daysince2 + foreign, data=panel,
               model="random", index = "id")
nwse <- coeftest(mod1for, vcov.=vcovNW)[,2]
pval <- coeftest(mod1for, vcov.=vcovNW)[,4]

screenreg(mod1for, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Posts per politician",
          custom.coef.map = list("invasion1"="Invasion",
                                 "woman1"="Woman",
                                 "invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "foreign1" = "Foreign",
                                 "daysince1"="Time",
                                 "daysince2" = "Time since invasion",
                                 "(Intercept)"="Intercept"),
          caption = "Effect of War on Politicians' Engagement, Controlling for Foreign Policy Elite",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:foreign",
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod1for$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469),
          reorder.gof = c(1, 3, 2))



# Table S4.1: Politician-Day 
polday1 <- polday[complete.cases(polday), ]
panel_sent_plm <- pdata.frame(polday1, index = c("id", "date"))
mod2zet <- plm(zet_m ~ daysince1 + invasion*woman + daysince2, data=panel_sent_plm,
               model = "random", index="id")
mod2nrc <- plm(nrc_m ~ daysince1 + invasion*woman + daysince2, data=panel_sent_plm,
               model="random", index="id")
mod2finn <- plm(afinn_m ~ daysince1 + invasion*woman + daysince2, data=panel_sent_plm,
                model="random", effect = "individual")
mod2bing <- plm(bing_m ~ daysince1 + invasion*woman + daysince2, data=panel_sent_plm,
                model="random", index="id")

nwse <- list(coeftest(mod2zet, vcov.=vcovNW)[,2],
             coeftest(mod2nrc, vcov.=vcovNW)[,2],
             coeftest(mod2finn, vcov.=vcovNW)[,2],
             coeftest(mod2bing, vcov.=vcovNW)[,2])
pval <- list(coeftest(mod2zet, vcov.=vcovNW)[,4],
             coeftest(mod2nrc, vcov.=vcovNW)[,4],
             coeftest(mod2finn, vcov.=vcovNW)[,4],
             coeftest(mod2bing, vcov.=vcovNW)[,4])

screenreg(list(mod2zet, mod2nrc, mod2finn, mod2bing),
          override.se = nwse,
          override.pvalues = pval,
          custom.header = list("Positive Sentiment" = 1:4),
          custom.model.names = c("SYUZHET", "NRC", "AFINN", "BING"),
          custom.coef.map = list("invasion1" = "Invasion",
                                 "woman1" = "Woman",
                                 "invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "daysince1" = "Time",
                                 "daysince2" = "Time since invasion",
                                 "(Intercept)" = "Intercept"),
          caption = "Effect of War on Positive Sentiment",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          label = "tab:sentiment",
          booktabs = T,
          stars=c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=c(round(sqrt(mod2zet$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod2nrc$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod2finn$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod2bing$ercomp[[1]][2]), digits=2)),
                                 "$N_{Politicians}$"=c(469, 469, 469, 469)),
          reorder.gof = c(1, 3, 2))


# Table S4.2: Effect of War on Topic Proportions (Binomial Logit)
topic_panel <- polday[complete.cases(polday), ]
modAidL <- glm(aid ~ daysince1 + invasion*woman + daysince2, data=topic_panel, weights = num.post,
               family = binomial(link="logit"))
modBoostL <- glm(boost ~ daysince1 + invasion*woman + daysince2, data=topic_panel, weights = num.post,
                 family=binomial(link="logit"))
modSecL <- glm(sec ~ daysince1 + invasion*woman + daysince2, data=topic_panel, weights = num.post,
               family = binomial(link="logit"))


a <- NeweyWest(modAidL, order.by= topic_panel$daysince1, verbose = T)
b <- NeweyWest(modBoostL, order.by= topic_panel$daysince1, verbose = T)
c <- NeweyWest(modSecL, order.by= topic_panel$daysince1, verbose = T)

nwse <- list(coeftest(modAidL, vcov=a)[,2], 
             coeftest(modBoostL, vcov=b)[,2], 
             coeftest(modSecL, vcov=c)[,2])
pval <- list(coeftest(modAidL, vcov=a)[,4], 
             coeftest(modBoostL, vcov=b)[,4], 
             coeftest(modSecL, vcov=c)[,4])

screenreg(list(modAidL, modBoostL, modSecL),
          override.se = nwse, override.pvalues = pval,
          custom.header = list("Topic proportion"=1:3),
          custom.model.names = c("Social aid",
                                 "Morale boost",
                                 "Security"),
          custom.coef.map = list("invasion1" = "Invasion",
                                 "woman1" = "Woman",
                                 "invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "daysince1" = "Time",
                                 "daysince2" = "Time since invasion",
                                 "(Intercept)" = "Intercept"),
          caption = "Effect of War on Topic Proportions (Binomial Logit)",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          include.aic=T,
          include.loglik=T,
          include.bic=F,
          booktabs = T,
          label = "tab:3topic_logit",
          booktabs = T,
          stars=c(0.1, 0.05, 0.01))


######################
##### Figure 4.1 #####
######################

post$label <- factor(post$label, levels = c("aid", "pride", "sec", "for", "leg", "rel", "rus", "other"))

plot_data <- post %>% ungroup() %>% select(woman, label, invasion) %>% 
  group_by(woman, invasion) %>%
  mutate(n.total = n()) %>% group_by(woman, invasion, label) %>%
  mutate(n = n(),
         prop = n/n.total,
         se = sqrt(prop*(1-prop)/n.total),
         ci.min = prop - 1.96*se,
         ci.max = prop + 1.96*se) %>% unique() %>% 
  mutate(across(c(prop, se), round, 5)) %>% group_by(woman, label) %>%
  mutate(diff = prop[invasion==1] - prop[invasion==0],
         diff.se = sqrt(se[invasion==1]^2 + se[invasion==0]^2))


s41 <- plot_data %>% filter(label %in% c("for", "leg", "rel", "rus", "other")) %>%
  ggplot(aes(x=label, y=diff, fill=woman)) + geom_bar(stat = "identity", position="dodge") +
  labs(x="Topic Labels", y="Change in Proportions", fill="Politician's \n Gender") + 
  scale_x_discrete(labels = c("Foreign affairs", "Legislation", "Religion", "Addressing Russia", "Other")) + 
  ylim(min = -.12, max = .1) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9"), labels=c("Men", "Women")) +
  geom_hline(yintercept=0, linetype=2, color="black") + 
  geom_errorbar(aes(ymin=diff-1.96*diff.se, ymax=diff+1.96*diff.se), width=.5, position=position_dodge(1)) +
  theme_classic() +
  theme(text = element_text(size = 15)) +
  theme(axis.text.x= element_text(size=10, angle=45, hjust=1))
s41

pdf("s41.pdf")
print(s41)
dev.off()





#### figure S4.2
## woman, before and after
s42women <- plot_data %>% filter(woman==1) %>%
  ggplot(aes(x=label, y=prop, fill=invasion)) + geom_bar(stat = "identity", position="dodge") +
  geom_errorbar(aes(ymin=ci.min, ymax=ci.max), width=.5, position=position_dodge(1)) +
  labs(x="Topic Labels", y="Proportion", fill="Invasion") + ggtitle("Women Politicians") +
  scale_x_discrete(labels = c("Social aid", "Morale boost", "Security", "Foreign affairs", "Legislation", "Religion", "Adressing Russia", "Other")) + 
  scale_fill_manual(values=c("#84c680", "#c280c6"), labels=c("Before", "After")) + theme_classic() +
  theme(text = element_text(size = 15)) +
  theme(axis.text.x= element_text(size=10, angle=45, hjust=1))
s42women

## man, before and after
s42men <- plot_data %>% filter(woman==0) %>%
  ggplot(aes(x=label, y=prop, fill=invasion)) + geom_bar(stat = "identity", position="dodge") +
  geom_errorbar(aes(ymin=ci.min, ymax=ci.max), width=.5, position=position_dodge(1)) +
  labs(x="Topic Labels", y="Proportion", fill="Invasion") + ggtitle("Men Politicians") +
  scale_x_discrete(labels = c("Social aid", "Morale boost", "Security", "Foreign affairs", "Legislation", "Religion", "Adressing Russia", "Other")) + 
  scale_fill_manual(values=c("#84c680", "#c280c6"), labels=c("Before", "After")) +  theme_classic() +
  theme(text = element_text(size = 15)) +
  theme(axis.text.x= element_text(size=10, angle=45, hjust=1)) 
s42men


pdf("s42women_prepost_full_ordered.pdf")
print(s42women)
dev.off()

pdf("s42men_prepost_full_ordered.pdf")
print(s42men)
dev.off()





# figure 5.1
plot_data <- polday %>% drop_na() %>%
  mutate(netreactions = react_m*num.post) %>%
  group_by(woman, date) %>% 
  summarise(date = date,
            n = sum(num.post),
            sum.react = sum(netreactions),
            m.react = sum.react/n,
            days = date - as.Date("2022-02-24")) %>% distinct()

s51 <- plot_data %>% 
  ggplot(aes(x = days, y = m.react, color=woman, shape=woman, linetype = woman)) +
  geom_point(size=2.5) + xlim(-100, 100) + 
  xlab("Days Since Invasion") + ylab("Reactions per Post") + 
  geom_vline(xintercept =0, linetype="dashed", color="red") +
  geom_smooth(data=subset(plot_data, days < 0), method="lm") + 
  geom_smooth(data=subset(plot_data, days > 0), method="lm") +
  theme_classic() +theme(text = element_text(size = 15)) + 
  scale_linetype_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c(3, 2)) +
  scale_color_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c("#E69F00", "#56B4E9")) +
  scale_shape_manual(name = "Politician's \nGender", labels=c("0"="Man", "1" = "Woman"), values=c(16, 17)) 

s51
pdf("s51react_full.pdf")
print(s51)
dev.off()




# Table S5.1: Effect of War on Public Reactions with Politician-Day Unit
panel_react_plm <- pdata.frame(polday, index = c("id", "date"))
mod3pd <- plm(react_m ~ daysince1 + invasion*woman + daysince2,
              model="random", index="id", data = panel_react_plm)
nwse <-coeftest(mod3pd, vcov.=vcovNW)[,2]
pval <-coeftest(mod3pd, vcov.=vcovNW)[,4]
screenreg(mod3pd, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Reactions",
          custom.coef.map =  list("invasion1"="Invasion",
                                  "woman1"="Woman",
                                  "invasion1:woman1" = "Invasion $\\times$ Woman",
                                  "daysince1"="Time",
                                  "daysince2" = "Time since invasion",
                                  "(Intercept)"="Intercept"),
          caption = "Effect of War on Public Reactions",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:reaction",
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod3pd$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469),
          reorder.gof = c(1, 3, 2))


#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Table S5.2: Effect of War on Public Reactions, Controlling for Office and Post
pfull <- pdata.frame(post, index = c("id", "dayseq"))
mod3tt <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2 + office + posttype,
              model="random", index="id", data = pfull)
nwse <-coeftest(mod3tt, vcov.=vcovNW)[,2]
pval <-coeftest(mod3tt, vcov.=vcovNW)[,4]
screenreg(mod3tt, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Reactions",
          custom.coef.map =  list("invasion1"="Invasion",
                                  "woman1"="Woman",
                                  "invasion1:woman1" = "Invasion $\\times$ Woman",
                                  "daysince1"="Time",
                                  "daysince2" = "Time since invasion",
                                  "(Intercept)"="Intercept"),
          caption = "Effect of War on Public Reactions",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod3tt$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469,
                                 "Office controls" = "\\checkmark",
                                 "Post type controls" = "\\checkmark"),
          reorder.gof = c(1, 3, 4, 5, 2),
          label = "table:twoway")


#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Table S5.3: Effect of War on Public Reactions for Different Reactions
mod3re <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2,
              model="random", index="id", data = pfull)
mod3like <- plm(Likes ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data = pfull)
mod3comm <- plm(Comments ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data = pfull)
mod3sh <- plm(Shares ~ daysince1 + invasion*woman + daysince2,
              model="random", index="id", data = pfull)
mod3love <- plm(Love ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data = pfull)
mod3an <- plm(Angry ~ daysince1 + invasion*woman + daysince2,
              model="random", index="id", data = pfull)

nwse <- list(coeftest(mod3re, vcov. = vcovNW)[,2],
             coeftest(mod3like, vcov. = vcovNW)[,2],
             coeftest(mod3comm, vcov. = vcovNW)[,2],
             coeftest(mod3sh, vcov. = vcovNW)[,2],
             coeftest(mod3love, vcov. = vcovNW)[,2],
             coeftest(mod3an, vcov. = vcovNW)[,2])
pval <- list(coeftest(mod3re, vcov. = vcovNW)[,4],
             coeftest(mod3like, vcov. = vcovNW)[,4],
             coeftest(mod3comm, vcov. = vcovNW)[,4],
             coeftest(mod3sh, vcov. = vcovNW)[,4],
             coeftest(mod3love, vcov. = vcovNW)[,4],
             coeftest(mod3an, vcov. = vcovNW)[,4])

screenreg(list(mod3re, mod3like, mod3comm, mod3sh, mod3love, mod3an),
          override.se = nwse, override.pvalues = pval,
          custom.model.names = c("Reactions", "Like", "Comment", "Share", "Love", "Angry"),
          custom.coef.map =  list("invasion1"="Invasion",
                                  "woman1"="Woman",
                                  "invasion1:woman1" = "Invasion $\\times$ Woman",
                                  "daysince1"="Time",
                                  "daysince2" = "Time since invasion",
                                  "(Intercept)"="Intercept"),
          caption = "Effect of War on Public Reactions for Different Reactions",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "table:types",
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=c(round(sqrt(mod3re$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod3like$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod3comm$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod3sh$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod3love$ercomp[[1]][2]), digits=2),
                                                                 round(sqrt(mod3an$ercomp[[1]][2]), digits=2)),
                                 "$N_{Politicians}$"=c(469, 469, 469, 469, 469, 469)),
          reorder.gof = c(1, 3, 2))



#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Table S5.4: Effect of War on Public Reactions with Topic Controls
pfull <- pfull %>% mutate(label = relevel(label, ref="other"))
mod3topic <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2 + label,
                 model="random", index="id", data = pfull)
nwse <-coeftest(mod3topic, vcov.=vcovNW)[,2]
pval <-coeftest(mod3topic, vcov.=vcovNW)[,4]
screenreg(mod3topic, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Reactions",
          custom.coef.map =  list("invasion1"="Invasion",
                                  "woman1"="Woman",
                                  "invasion1:woman1" = "Invasion $\\times$ Woman",
                                  "daysince1"="Time",
                                  "daysince2" = "Time since invasion",
                                  "labelaid" = "Social aid",
                                  "labelpride" = "Morale boost",
                                  "labelsec" = "Security",
                                  "labelfor" = "Foreign affairs",
                                  "labelleg" = "Legislation",
                                  "labelrel" = "Religion",
                                  "labelrus" = "Addressing Russia",
                                  "(Intercept)"="Intercept"),
          caption = "Effect of War on Public Reactions with Topic Controls",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod3topic$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469),
          reorder.gof = c(1, 3, 2),
          label = "table:topic")


#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Table S5.5: Effect of War on Public Reactions with Region Controls
mod3reg <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2 + region,
               model="random", index="id", data = pfull)
nwse <-coeftest(mod3reg, vcov.=vcovNW)[,2]
pval <-coeftest(mod3reg, vcov.=vcovNW)[,4]
screenreg(mod3reg, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Reactions",
          custom.coef.map =  list("invasion1"="Invasion",
                                  "woman1"="Woman",
                                  "invasion1:woman1" = "Invasion $\\times$ Woman",
                                  "daysince1"="Time",
                                  "daysince2" = "Time since invasion",
                                  "(Intercept)"="Intercept"),
          caption = "Effect of War on Public Reactions with Region Controls",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod3reg$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=469,
                                 "Region controls" = "\\checkmark"),
          reorder.gof = c(1, 3, 4, 2),
          label = "table:reg")


#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Table S5.6: Effect of War on Public Reactions, Controlling for Followers
pfull1 <- pfull %>% filter(nfollower != "N/A")
pfull1$nfollower <- as.numeric(pfull1$nfollower)
mod3fol <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2 + nfollower,
               model="random", index="id", data = pfull1)
nwse <-coeftest(mod3fol, vcov.=vcovNW)[,2]
pval <-coeftest(mod3fol, vcov.=vcovNW)[,4]
screenreg(mod3fol, override.se = nwse, override.pvalues = pval,
          custom.model.names = "Reactions",
          custom.coef.map =  list("invasion1"="Invasion",
                                  "woman1"="Woman",
                                  "invasion1:woman1" = "Invasion $\\times$ Woman",
                                  "daysince1"="Time",
                                  "daysince2" = "Time since invasion",
                                  "nfollower" = "Follower",                                  
                                  "(Intercept)"="Intercept"),
          caption = "Effect of War on Public Reactions, Controlling for Followers",
          caption.above = T,
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$\\hat{\\sigma}_\\text{id}$"=round(sqrt(mod3fol$ercomp[[1]][2]), digits=2),
                                 "$N_{Politicians}$"=274),
          reorder.gof = c(1, 3, 2),
          label = "tab:follower")



#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Table S5.7: Popular Posts by Topic
post$label <- factor(post$label, levels = c("aid", "pride", "sec", "for", "leg", "rel", "rus", "other"))
full_man <- post[post$woman==0, ]
full_woman <- post[post$woman==1, ]
top100int_man <- full_man[c(order(full_man$totalreactions, decreasing = T)[1:100]),]
top100int_woman <- full_woman[c(order(full_woman$totalreactions, decreasing = T)[1:100]),]
xtabs(~ label, data=top100int_woman)
xtabs(~ label, data=top100int_man)

# Table S5.8: Average Sentiment Scores of Popular Posts
round(colMeans(top100int_woman[, c("zet", "nrc", "bing", "afinn")]), digits=2)
round(colMeans(top100int_man[, c("zet", "nrc", "bing", "afinn")]), digits=2)




#######################
###### Placebos #######
#######################

# Table S6.1: Placebo Invasions and Posting Behavior
prepanel <- polday %>% filter(invasion ==0)

# 12-28
panel1228 <- prepanel %>% 
  mutate(invasion = factor(ifelse(date < as.Date("2021-12-28"), 0, 1)),
         daysince2 = ifelse(date - as.Date("2021-12-27") < 0, 0, date - as.Date("2021-12-27"))) %>%
  group_by(id) %>% mutate(daysince1 = seq.int(n()))
prepanel1 <- pdata.frame(panel1228, index = c("id", "date"))
placebo1 <- plm(num.post ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data=prepanel1)

# 02-22
panel222 <- prepanel %>% 
  mutate(invasion = factor(ifelse(date < as.Date("2022-02-22"), 0, 1)),
         daysince2 = ifelse(date - as.Date("2022-02-21") < 0, 0, date - as.Date("2022-02-21"))) %>%
  group_by(id) %>% mutate(daysince1 = seq.int(n()))
prepanel1 <- pdata.frame(panel222, index = c("id", "date"))
placebo2 <- plm(num.post ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data=prepanel1)

# 02-17
panel217 <- prepanel %>% 
  mutate(invasion = factor(ifelse(date < as.Date("2022-02-17"), 0, 1)),
         daysince2 = ifelse(date - as.Date("2022-02-16") < 0, 0, date - as.Date("2022-02-16"))) %>%
  group_by(id) %>% mutate(daysince1 = seq.int(n()))
prepanel1 <- pdata.frame(panel217, index = c("id", "date"))
placebo3 <- plm(num.post ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data=prepanel1)

# 02-10
panel210 <- prepanel %>% 
  mutate(invasion = factor(ifelse(date < as.Date("2022-02-10"), 0, 1)),
         daysince2 = ifelse(date - as.Date("2022-02-09") < 0, 0, date - as.Date("2022-02-09"))) %>%
  group_by(id) %>% mutate(daysince1 = seq.int(n()))
prepanel1 <- pdata.frame(panel210, index = c("id", "date"))
placebo4 <- plm(num.post ~ daysince1 + invasion*woman + daysince2,
                model="random", index="id", data=prepanel1)

nwse <- list(coeftest(placebo1, vcovNW)[,2],       
             coeftest(placebo2, vcovNW)[,2],
             coeftest(placebo3, vcovNW)[,2],
             coeftest(placebo4, vcovNW)[,2])
pval <- list(coeftest(placebo1, vcovNW)[,4],       
             coeftest(placebo2, vcovNW)[,4],
             coeftest(placebo3, vcovNW)[,4],
             coeftest(placebo4, vcovNW)[,4])

screenreg(list(placebo1, placebo2, placebo3, placebo4),
          override.se = nwse, override.pvalues = pval,
          custom.header = list("Posts per politician" = 1:4),
          custom.model.names = c("2021-12-28", "2022-02-22", "2022-02-17", "2022-02-10"),
          custom.coef.map = list("invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "(Intercept)"="Intercept"),
          caption = "Different Placebo Invasions on Gendered Effect of War",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=T,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:placebo1",
          stars = c(0.1, 0.05, 0.01),
          custom.gof.names = c('$N_{Observations}$'),
          custom.gof.rows = list("$N_{Politicians}$"=c(469, 469, 469, 469)),
          reorder.gof = c(2,1))

# Table S6.2: Placebo Invasions and Other Outcomes
# filter to pre-invasion posts
# get the running variable value for the invasion date
unique(polday$daysince1[polday$date=="2022-02-24"]) # 116
prefull <- post %>% filter(daysince1 < 116)
prefull <- prefull %>% group_by(id, daysince1) %>%
  mutate(seq = seq.int(n()),
         dayseq = paste0(daysince1,".",seq))

# Sentiments
# 2021-12-28
unique(polday$daysince1[polday$date=="2021-12-28"]) # 58
unique(polday$daysince1[polday$date=="2021-12-27"]) # 57
prefull1 <- prefull %>% 
  mutate(invasion = as.factor(ifelse(daysince1 < 58, 0, 1)),
         daysince2 = ifelse(daysince1 -57<0, 0, daysince1-57))
pfull1 <- pdata.frame(prefull1, index = c("id", "dayseq"))
placebo1.1 <- plm(zet ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull1)
placebo1.2 <- plm(nrc ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull1)
placebo1.3 <- plm(afinn ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull1)
placebo1.4 <- plm(bing ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull1)

# 2.22
unique(polday$daysince1[polday$date=="2022-02-22"]) # 114
unique(polday$daysince1[polday$date=="2022-02-21"]) # 113
prefull2 <- prefull %>% 
  mutate(invasion = as.factor(ifelse(daysince1 < 114, 0, 1)),
         daysince2 = ifelse(daysince1 - 113<0, 0, daysince1-113))
pfull2 <- pdata.frame(prefull2, index = c("id", "dayseq"))
placebo2.1 <- plm(zet ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull2)
placebo2.2 <- plm(nrc ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull2)
placebo2.3 <- plm(afinn ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull2)
placebo2.4 <- plm(bing ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull2)

# 2.17
unique(polday$daysince1[polday$date=="2022-02-17"]) # 109
unique(polday$daysince1[polday$date=="2022-02-16"]) # 108
prefull3 <- prefull %>% 
  mutate(invasion = as.factor(ifelse(daysince1 < 109, 0, 1)),
         daysince2 = ifelse(daysince1 - 108<0, 0, daysince1-108))
pfull3 <- pdata.frame(prefull3, index = c("id", "dayseq"))
placebo3.1 <- plm(zet ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull3)
placebo3.2 <- plm(nrc ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull3)
placebo3.3 <- plm(afinn ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull3)
placebo3.4 <- plm(bing ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull3)

# 2.10
unique(polday$daysince1[polday$date=="2022-02-10"]) # 102
unique(polday$daysince1[polday$date=="2022-02-09"]) # 101
prefull4 <- prefull %>% 
  mutate(invasion = as.factor(ifelse(daysince1 < 102, 0, 1)),
         daysince2 = ifelse(daysince1 - 101<0, 0, daysince1-101))
pfull4 <- pdata.frame(prefull4, index = c("id", "dayseq"))
placebo4.1 <- plm(zet ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull4)
placebo4.2 <- plm(nrc ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull4)
placebo4.3 <- plm(afinn ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull4)
placebo4.4 <- plm(bing ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull4)

# Topics
# 12.28
topic1 <- panel1228[complete.cases(panel1228), ]
topic1 <- pdata.frame(topic1, index = c("id", "date"))
placebo1.5 <- plm(aid ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic1)
placebo1.6 <- plm(boost ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic1)
placebo1.7 <- plm(sec ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic1)

# 2.22 
topic2 <- panel222[complete.cases(panel222), ]
topic2 <- pdata.frame(topic2, index = c("id", "date"))
placebo2.5 <- plm(aid ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic2)
placebo2.6 <- plm(boost ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic2)
placebo2.7 <- plm(sec ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic2)

# 2.17
topic3 <- panel217[complete.cases(panel217), ]
topic3 <- pdata.frame(topic3, index = c("id", "date"))
placebo3.5 <- plm(aid ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic3)
placebo3.6 <- plm(boost ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic3)
placebo3.7 <- plm(sec ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic3)


# 2.10 
topic4 <- panel210[complete.cases(panel210), ]
topic4 <- pdata.frame(topic4, index = c("id", "date"))
placebo4.5 <- plm(aid ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic4)
placebo4.6 <- plm(boost ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic4)
placebo4.7 <- plm(sec ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=topic4)


#############################################
##### REDACTED DUE TO DATA RESTRICTIONS #####
#############################################
# Reactions
placebo1.8 <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull1)
placebo2.8 <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull2)
placebo3.8 <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull3)
placebo4.8 <- plm(totalreactions ~ daysince1 + invasion*woman + daysince2,
                  model="random", index="id", data=prefull4)


screenreg(list(coeftest(placebo1.1, vcovNW),
               coeftest(placebo1.2, vcovNW),
               coeftest(placebo1.3, vcovNW),
               coeftest(placebo1.4, vcovNW),
               coeftest(placebo1.6, vcovNW),
               coeftest(placebo1.7, vcovNW),
               coeftest(placebo1.5, vcovNW)),
               #coeftest(placebo1.8, vcovNW)),
          custom.header = list("Sentiment" = 1:4,
                               "Topics" = 5:7),
                               #"Reactions" = 8),
          #custom.model.names = c("2021-12-28", "1 Day Before", "1 Week Before", "2 Weeks Before"),
          custom.model.names = c("ZET", "NRC", "BING", "AFINN", "Social aid", "Morale boost", "Security"),
          custom.coef.map = list("invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "(Intercept)"="Intercept"),
          caption = "Different Placebo Invasions on Gendered Effect of War for Different DVs",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=F,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:placebo2",
          scalebox = 0.8,
          stars = c(0.1, 0.05, 0.01))

screenreg(list(coeftest(placebo2.1, vcovNW),
               coeftest(placebo2.2, vcovNW),
               coeftest(placebo2.3, vcovNW),
               coeftest(placebo2.4, vcovNW),
               coeftest(placebo2.6, vcovNW),
               coeftest(placebo2.7, vcovNW),
               coeftest(placebo2.5, vcovNW)),
               #coeftest(placebo2.8, vcovNW)),
          custom.header = list("Sentiment" = 1:4,
                               "Topics" = 5:7),
                               #"Reactions" = 8),
          #custom.model.names = c("2021-12-28", "1 Day Before", "1 Week Before", "2 Weeks Before"),
          custom.model.names = c("ZET", "NRC", "BING", "AFINN", "Social aid", "Morale boost", "Security"),
          custom.coef.map = list("invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "(Intercept)"="Intercept"),
          caption = "Different Placebo Invasions on Gendered Effect of War for Different DVs",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=F,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:placebo2",
          scalebox = 0.8,
          stars = c(0.1, 0.05, 0.01))

screenreg(list(coeftest(placebo3.1, vcovNW),
               coeftest(placebo3.2, vcovNW),
               coeftest(placebo3.3, vcovNW),
               coeftest(placebo3.4, vcovNW),
               coeftest(placebo3.6, vcovNW),
               coeftest(placebo3.7, vcovNW),
               coeftest(placebo3.5, vcovNW)),
               #coeftest(placebo3.8, vcovNW)),
          custom.header = list("Sentiment" = 1:4,
                               "Topics" = 5:7),
                               #"Reactions" = 8),
          custom.model.names = c("ZET", "NRC", "BING", "AFINN", "Social aid", "Morale boost", "Security"),
          custom.coef.map = list("invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "(Intercept)"="Intercept"),
          caption = "Different Placebo Invasions on Gendered Effect of War for Different DVs",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=F,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:placebo2",
          scalebox = 0.8,
          stars = c(0.1, 0.05, 0.01))

screenreg(list(coeftest(placebo4.1, vcovNW),
               coeftest(placebo4.2, vcovNW),
               coeftest(placebo4.3, vcovNW),
               coeftest(placebo4.4, vcovNW),
               coeftest(placebo4.6, vcovNW),
               coeftest(placebo4.7, vcovNW),
               coeftest(placebo4.5, vcovNW)),
               #coeftest(placebo4.8, vcovNW)),
          custom.header = list("Sentiment" = 1:4,
                               "Topics" = 5:7),
                               #"Reactions" = 8),
          custom.model.names = c("ZET", "NRC", "BING", "AFINN", "Social aid", "Morale boost", "Security"),
          custom.coef.map = list("invasion1:woman1" = "Invasion $\\times$ Woman",
                                 "(Intercept)"="Intercept"),
          caption = "Different Placebo Invasions on Gendered Effect of War for Different DVs",
          caption.above = T,
          include.rsquared = F,
          include.adjrs = F,
          include.nobs=F,
          include.rs=F,
          include.variance=F,
          include.deviance=F,
          booktabs = T,
          label = "tab:placebo2",
          scalebox = 0.8,
          stars = c(0.1, 0.05, 0.01))
